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Abstract 

We describe the implementation of the shearingbox approach into the Godunov- 
type central-upwind/constraint-transport magnetohydrodynamics code NIRVANA. 
This will allow for applications which require sheared-periodic boundary conditions 
as typically used in local Cartesian simulations of differentially rotating systems. We 
present the algorithm in detail and discuss necessary modifications in the numerical 
fluxes in order to preserve conserved quantities and to fulfill other analytical con- 
straints as good as seem feasible within the numerical scheme. We check the source 
terms which come with the shearingbox formulation by investigating the conserva- 
tion of the epicyclic mode energy. We also perform more realistic simulations of the 
magneto-rotational instability with initial zero-net-flux vertical magnetic field and 
compare the obtained stresses and energetics with previous non-conservative results 
exploring the same parameter regime. 
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1 Introduction 

Simulations aiming to study the local behavior of (magneto-) hydrodynamic 
flows usually invoke periodic boundary conditions. This accounts for the un- 
boundedness of the flow and prevents the system from obtaining information 
about its surface. However, in the case of a shearing background flow period- 
icity is not strictly applicable. Typical scenarios for this include differentially 
rotating gaseous accretion disks, which are encountered in many astrophysical 
contexts. Such systems (if sufficiently ionized) are found to be unstable to 
the so-called magneto-rotational instability (MRI) and have been successfully 
modeled during the past decade (see [2f3] and references therein). Consider- 
able progress in the field has been made using the local shearingbox approach, 
a detailed description of which can be found e.g. in pp. 
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The two key issues to be addressed in the design of shearingbox boundary 
conditions are (1) the time-dependent shifted-periodic mapping of the depen- 
dent variables of mass density, momentum, energy and, in case of magneto- 
hydrodynamic flows, the magnetic field and (2) the offset introduced in the 
velocity component parallel to the shearing background. The latter represents 
the global shear across the domain which is forced constant in time. Both 
of these features and the numerical treatment of the resulting source terms 
arising in the local approximation have profound implications on the accuracy 
properties of the shearingbox advection/induction scheme. 

In this paper we will deal with these implications in detail. First, section 2 gives 
a compact recollection of basic facts about the MHD scheme in NIRVANA 
necessary for a understanding of the following shearingbox implementation 
which is then described in Section 3. Section 4 presents test cases to validate 
our approach. 



2 Review of the MHD scheme 



Although the NIRVANA code can account for dissipative effects (viscosity, 
diffusion, thermal heat conduction), multi-scale phenomena (using adaptive 
mesh refinement, see Ziegler [T5]) as well as self-gravity (multigrid-type Pois- 
son solver, see Ziegler [13]), we want to focus here on ideal MHD in three space 
dimensions described by the following set of coupled equations 



<9u + dP + dP + dP g (u) 

dt dx dy dz 



dB 

~dt 



+ VxE = 0. 



(1) 
(2) 



where u = (g,m x ,m y ,m z ,e) T . The source term S arises from non-inertial 
forces in the local frame (rotating with O relative to the inertial frame) and 
will be discussed in detail in section 13.31 The flux functions P, P. P are 
defined as: 
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We adopt a notation where the magnetic permeability is set to unity. The 
equations are to be supplemented by the zero- divergence constraint V-B = 
for the magnetic field B and an equation of state e = e + \qv 2 + |B 2 with 
e = p/{l ~ !)• Other symbols have their usual meaning with g denoting the 
mass density, e the total energy density, e the thermal energy density, p the 
pressure, v the velocity, m = gv the momentum, E = — v x B the electric 
field and 7 the ratio of specific heats. 

The numerical method used to solve the above system is a hybrid scheme com- 
bining the second-order version of the Godunov-type central-upwind scheme 
of Kurganov et al. [8] with the constraint transport (CT) technique for diver- 
gence free magnetic field evolution (see e.g. Evans & Hawley [5], Toth [10J). 
The original description of this hybrid approach can be found in Ziegler 
Here, we compactly summarize the relevant knowledge necessary for construct- 
ing the shearingbox-implementation. 

Let the computational domain be subdivided in uniform Cartesian cells Cij t k = 
[x. i,x. 1] x [y . i,y, 1] x [z i,z 1] with center (xi,yj,Zk) and spacings Sx, 

t- 2 i+ 2 j- 2 j+ 2 k- 2 k+ 2 

Sy and 8z. Then, the scheme for Eq. ([T]) in semi-discrete form reads 



F»\ (t)-Pf! .it) F? , (f)-F? (f) 
di UljA) = Tx Ty 



F* (t)-F* ,{t) 
Jz 



+ S(ui J; k (t)) 



(6) 



where the cell-centered Ujj^ is an approximation to the volume-averaged state 
u in cell Ci^^ and the numerical fluxes are given by: 
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with w= (u, B) T . We apply the third-order strong-stability-preserving Runge- 
Kutta scheme in [7j to integrate Eq. (jBJ) in time (see also Eq. ( |T9l) ). At any 
stage in the Runge-Kutta scheme the evaluation of the flux functions at cell 
faces requires point values for the variables left / right to a cell interface which 
are obtained by piecewise linear reconstruction from Wjj^ applying the TVD 
slope-limiter of van Leer [9]. These point values are characterized in the fol- 
lowing by superscripts W at cell interface x = x. i, Eat x = x. i, S at 

y — y. i, N at y — y. i, B at z — z i and T at z — z i, respectively With 

j~2 j+2 «- 2 K+2 

such notation w^ ljjfc stands for the reconstructed value within cell C i+ ij^ 
at x = x i and wr,- & for the value within cell ,• t at the same position. 

In contrast to the cell- centered Uij,k we use staggering for the magnetic field 
i.e. the face- centered finite-volume approximations to the face-averaged com- 
ponents are B , i , B , i and i, respectively. It follows that in 

the reconstruction of B at E,W position the x-component is already properly 
located, i.e. (B x )ff% — B i . , whereas other components are first recon- 

structed in x-direction and then spatially averaged in the transverse direction. 
Reconstruction of B at the other locations N,S,B,T is done in analogous fash- 
ion. 



The quantities a 1 * 1 in the flux formulae define the maximum (plus sign) respec- 
tive minimum (minus sign) wave-propagation-direction-sensitive speed at the 
interface x 



if 2 ' 



a+ =m&x{(v x + Cf)JJ 1J)fc , (v x + Cf)£ >fe , 0} , 
a" =min{(v x - ct)^. 1Jtk , (v x - c f )g ifc , 0} , 



where Cf = (c| + c^) 1//2 is an upper limit for the fast magnetosonic speed 
including the sound speed c s = (7p/f?) 1//2 and Alfven-speed ca = (B 2 /^) 1 / 2 . 
The corresponding speeds in the ?/(z)-direction are denoted by b^ 1 (c ). 

The above method can be formally applied to the induction equation (j2J). 
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Writing the curl of the electric field as divergence of an antisymmetric tensor 



-V x E = V- 



( E z -E v \ 
—E z E x 

\ By —E x 



and utilizing the abbreviations e x = (0, — E z , E y ) T , e y = (E z , 0, — E X ) T , and 
e z = (—E y ,E x ,0) T , in analogy to Eqs. (7)-(9) one can derive electric field 
fluxes 
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These are, like the F-fluxes, defined at cell faces. As before, at any stage 
in the Runge-Kutta scheme, the electric field function e x (e v , e z ) have to 
be evaluated at E,W (N,S,T,B) cell interfaces using reconstructed data, e.g., 
^(wfy = {-E,)f Jtk = -(v x )f Jtk (B v )f Jtk • (r^JiUljj, Our CT scheme 
for the staggered magnetic field components finally reads in semi-discrete form 



a . . i . i , — E i . i E i i— E i i 

11° i- 1 ■ u = r 1 r ' l 1 ^ 



di 
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(15) 



where are edge-centered approximations to 

x \ l ->3~2'' k ~2 y\ l ~ 2>3> 2 z \ % ~2^~2^ 

the edge-averaged electric field components. These are computed by compo- 
sition of the electric field fluxes obtained from the Godunov-central scheme, 
i.e., 



E, i i = -(-G y 1 -G y 1 +G Z i+G z i), (16) 

x\i,j-2,k-2 4 V z\i,j~2,k z\i,j-^,k-l y\i,3,k-^ y\i,j-l,k-^ J 
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?/|«— 2 >J , ': fc — 2 4 ^[«-2'i' fc zK-^J.^^ 1 x\i,j,k-^ x\i— l,j,k—^ 

£,1.1 =7(-G r i - i +G y J +G y ! ). (18) 

z|j-2>-?-2' fc 4 y|i-2,i,fc tf|*-2,J-l,fe x\i,j-^,k x\i-l,j-^,k' 

It is now easy to show from f|T3l) - f|T5|) that 



x\i{--^,j,k x\t-^J,k 

5x 



52 



0. 



Thus, if (V-B)jj'fc = initially, the system evolves divergence-free since this 
condition is fulfilled after each stage in the multi-stage Runge-Kutta scheme. 

The full system of ODEs (jSJ), (IT3"|) - (fT5|) including the source term is solved with 
the strong-stability-preserving third-order Runge-Kutta integration scheme 
given by (dropping cell indices and writing the rhs of Eq. O as ~Lp + S and 
(HHD (USD compactely as |B = h E ) 



u« =u n + 5t{L n F + S n ) 
B (1) =B"' + 5tL^, 

u( 2 ) = ^ + V) + Ift(L«+SW) 
4 4 4 v * ' 

B (2) = -B" + 1b (1) + iftLgJ , 
4 4 4 s ' 

_ rt+1= l n+ 2 (2) + 2 t(L ( 2) + g(2)) 

o o o 

gn+1 = lgn 2 ( 2) 2 g) 

3 3 3 
with the time-step 5i = t n+1 — t n restricted by the CFL condition. 



3 Shearingbox implementation 



Let us now consider a magnetohydrodynamic flow with a background profile 
of the form v = (0, v y = sx, 0), with (linear) shear-parameter s. The compu- 
tational domain is a Cartesian box with dimensions L x , L y , L z and periodic 
boundaries in y— and z— direction. To account for the background shear we 



6 




Fig. 1. Mapping of ghost zone values for a 2D shear ingsheet. Reconstruction is 
indicated by highlighted cells. Arrows illustrate the need for a global velocity offset 
to match the velocity profile of the neighboring domain. 

introduce shearingbox boundary conditions in x— direction. These can be ex- 
pressed mathematically in the forrrPI (e.g. pQ): 



f{x,y,z) i-> f(x±L x ,yTwt,z), f E {g, m x , m z , e} (20) 

m y (x, y, z) i-> m y (x ±L x ,y^p wt, z) =f gw, (21) 

where w = —sL x represents the global velocity-offset across the box, as can be 
seen in Figured! Because NIRVANA being a conservative code evolves the total 
energy (rather than the internal energy e) we supplement the corresponding 
relation for the total energy density 



e(x, y, z) \— y e(x ± L x , y =p wt, z) =p m y w + \gw s 



[22) 



which can be easily derived by separating the rriy/(2g) part from the total 
energy and then using (I2T!) . 

As the y-coordinate of above mappings varies continously in time, there is 
some kind of interpolation necessary to map ghostzone values on a finite grid 
(see e.g. pQ for details). For our implementation, the same piecewise linear 
reconstruction is used as for the numerical scheme. 



Due to the shifted periodicity and the additional interpolation, ghost cells on 
opposing sides of the domain along with their adjoint regular zones are not re- 



1 Position arguments are being suppressed for extra variables in eqs. ([2T]) . ([22]) . 
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dundant in the same sense as for strictly periodic boundary conditions. Let us 
elucidate on that: In the case of regular periodicity every pair of cell and ghost 
cell (separated by the domain boundary) has an identical counterpart at the 
opposite boundary. Therefore the total flux across the interface is conserved 
to machine accuracy. 

For shifted periodicity at a given time t this does no longer hold. Thus the ad- 
ditional truncation error connected to the interpolation can lead to the buildup 
of considerable deviations in conserved quantities, as is shown in section 14.11 
In the following we want to suggest a method to correct for this. Since strict 
conservation only applies as long as there is no source on the rhs of (jTJ, in 
the local shearingbox approach with source term (1331) / (1341) only the total z- 
momentum and mass are exactly conserved. We note that it is nevertheless 
desirable to preserve all variables as good as possible under advection. 

3. 1 Conservation of hydrodynamic variables 

As the numerical fluxes are nonlinear functions of the conserved quantities 
any form of interpolation for the ghost cells will lead to some small inconsis- 
tency in the fluxes. This can be avoided by matching the computed x-fluxes 
at the sheared domain boundaries. It is straightforward to map fluxes not 
containing m y . The quantities related to those fluxes, i.e. mass density, x- 
and ^-momentum are then conserved to machine precision with respect to 
advection. For the ^-momentum flux and total energy flux, applying the map- 
pings (|2"Tj) / (J22J) to the third and fifth component of the flux function Q yields: 



in nice analogy to the original mappings. We want to point out that these rela- 
tions consistently include the magnetic part of the fluxes due to the Lorentz- 
force, although this is not directly visible in the above notation. Also, the 
modifications to the fluxes can be solely expressed in terms of fluxes and the 
velocity offset w. 

The recipe can then directly be carried over to the numerical fluxes ([7]): while 
the first term is just a linear combination of the flux-function (evaluated at 
the adequate positions) one has to apply (l2~il / (1221) once again for the second 
term, which then gives: 



f X ( m y) >-> f X ( m y) T f X {0) w j 
f(e)^f(e)T/» + if(^ 2 , 



(23) 
(24) 



K) =F*i K) T F£i (q)w, 



(25) 
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F :±^ e ) = p i±~ 3h ^ =f ^ JJfc K) w + ^j g)w2 ' (26) 



where the hat stands for the piecewise linear interpolation procedure used, 
and tilde marks the corresponding indices of the zones to map from. 

Relation ( 1241 expresses the fact that the total energy within the shearingbox 
is not conserved but can be altered by angular momentum transport through 
the radial boundaries. In integral formulation this corresponds to Eq. (8) (cast 
into our notation) of Hawley et al. pQ: 

dT r 

— = w J dy dz (gvjvy - B x B y ) (27) 
ax 



where V is the volume integral of total energy and 5v y = v y + sx expresses the 
perturbed velocity!" 2 ! When correcting the energy- and momentum-fluxes ac- 
cording to Eqs. (|23|) . (1241) one can satisfy this property more accurately which 
allows to trace the detailed evolution of energetics, e.g. in MRI-simulations. 



3.2 Conservation of magnetic flux 



Applying GauB' theorem to the integral form of the induction equation, one 
can show that the azimuthal field (due to the shear) grows linearly with the 
net radial magnetic field-flux through the radial boundaries: 

B -g>=-$9j*y*B,. (28) 

dX 



For zero net radial field the mean magnetic flux through the shearingbox is 
conserved. Our implementation satisfies this condition to machine accuracy 
for the x- and z-component of the magnetic field and to truncation error for 
the y-component (see section H~T]) . This corresponds to the amount of precision 
reported by [TH6] . 

The seminal paper of Hawley et. al [lj discusses this topic rather briefly, men- 
tioning that mapping the electromotive forces at the sheared interfaces con- 
serves the vertical field to roundoff error. Furthermore, the authors state that 
the spurious vertical field, that arises otherwise, is negligible as the associ- 
ated MRI growth rates are unresolved. Considering the differences between 
MRI-simulations with and without zero net vertical field (see e.g. [6]) leaves 

2 Expanding 5v y yields the additional w 2 term in 
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this somewhat questionable in the sense that a spurious field might lead to 
an overestimation of magnetic stresses in the latter case. One might further 
investigate this by a series of MRI-simulations with initial vertical field of the 
form B z = B sm(2TTx/L x ) + Bi, where \Bi\ <C \Bq\, i.e., with a controlled 
"spurious" field. Although from the physical point of view an exactly vanish- 
ing vertical flux is rather a question of academic nature, the issue of enhanced 
turbulent transport for net vertical field is still standing. 

In our implementation we apply additional boundary conditions to the electric 
field fluxes G to assure conservation of mean magnetic fields. Due to the 
different staggering of magnetic field components, which are face-centered, the 
method has to be applied for the fluxes in all three directions. The velocity 
offset w enters the fluxes via the electromotive force: 

E(x, y, z) h-> E(x ± L x ,y=pwt,z) ±w y xB (29) 



For the x(^)-direction only the x(z)-component of the magnetic field is needed 
to evaluate the additional second term. As the staggering of the normal field 
component coincides with that of the fluxes no reconstruction is needed in 
this case. Employing B^ i jk = S^+ij,* (respectively B^. jk = B% iJjk+1 ) the 
characteristic velocities cancel out and we yield 



G% . = Gfi _ =f w yB i 



G z i 



G 2 



i TwyB 



z\i,jj,k+: 



(30) 
(31) 



for the numerical electric fluxes in x/z-direction. For the remaining direction 
matters are a bit more complicated and we need the full reconstruction to be 
consistent with the underlying numerical scheme: 

G y ! = G y ! ±w — : — (32) 



To conclude the discussion of the boundary conditions we want to remark that 
the reconstruction of the magnetic field at the sheared interfaces preserves the 
V-B = constraint to roundoff error (see section Ei~T]) . The actual code is pub- 
licly available on the internet at www.aip.de/~gressel/index.php?id=code 
and can be embedded into the original NIRVANA package. 

The current implementation supports distributed memory parallelism in the 
form of blockwise domain decomposition along the z-coordinate. In this case all 
the information needed to reconstruct the boundary values is available locally, 
i.e., without MPI communications. Additional block distribution along the x- 



10 



and y-coordinates would in principle improve the surface to volume ratio but 
would also add communication overhead. This is discussed in detail in sections 
4.3 and 4.4 of [4J. Their Figure 5 shows that block distribution along the 
x-coordinate is inefficient and performance gains for additional distribution 
along the y-direction are moderate. They argue that the poor performance for 
x-distribution might be a cache issue with the Fortran column major ordering, 
which leaves this result marginally meaningful for our implementation. In 
the case of vertically elongated boxes, like they are common for simulations 
in galactic environments, distribution along z even outperforms y- and y-z- 
distributions. 



3.3 Source terms 

The second substantial ingredient of the shearingbox formalism are the source 
terms in the so called Hill system. This approximation is based on the local 
expansion of the equations of motion resulting in a tidal force 2qQ 2 x where 
q = dlnfi/dlni? represents the (angular) shear parameter for a differential 
rotation of the form Q(R) oc R~ q . Together with the Coriolis force the source 
terms for the momentum- and energy equation are: 

S(m) = -2gilz x v + 2gqil 2 xk = -2oQi x (v + qttxy) (33) 
5(e) = +2ett 2 qxx ■ v. (34) 

Equation (1331 shows that the momentum source terms can be combined and 
act as an effective^ Coriolis force on the perturbed velocity 5v y = v y + qVLx. 
This would in principle allow for an exact Coriolis-update in the form of 
an analytic rotation. Numerical tests, however, show that such an update 
is not suitable for the multi-stage Runge-Kutta integration scheme. Because 
operator-splitting is not favorable (for reasons described in section 14.21 below) 
we decide to implement the source terms unsplit, i.e., as explicit forces within 
the Runge-Kutta time integration scheme (see eq. (fT51) ). 

Gardiner & Stone [5] stress the importance of conserving the energy contained 
in the epicyclic mode. This ideally conserved quantity can be derived from the 
energy budget in the limit of inviscid flow. The general expression reads: 

I 2 ( ^ ) 1 2 1 i 2 ^ 2 i 



3 In contrast to the description in the Lagrangian frame of reference the angular 
velocity = is not a function of radius here. 
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resolution 




24 2 


32 2 


48 2 


64 2 


96 2 


mass 


a 


1.3371e-07 


1.2407e-07 


6.0259e-08 


1.2674e-08 


8.5488e-09 




b 


4.8018e-16 


1.1181e-15 


8.0027e-16 


1.4405e-15 


1.6005e-15 



Table 1 



Convergence study of the relative error in total mass: (a) without special treatment, 
(b) with the method described above. 

with k the epicyclic frequency. This formally looks like a kinetic energy but 
also includes the potential energy with respect to the epicyclic displacement. 
We have found that it is important to implement the source terms in an 
unsplit fashion to avoid oscillations in this energy that would otherwise arise 
from systematic splitting errors. This is discussed in more detail in section H~2l 
below. 



4 Test cases 

We have validated our implementation with various simple advection tests, i.e., 
with source terms switched off, and have checked the mentioned conservation 
properties. 

4-1 Advection tests 

As an example Table [1] shows the maximum relative error in the total mass 
contained in our simulation box. Without special treatment of the fluxes we 
find an error that decreases with increasing resolution, i.e., has the form of a 
truncation error. With our modified treatment the total mass is conserved to 
roundoff error. Albeit not explicitely shown in Table [IJ this is also true for the 
momentum and total energy. 

As discussed in section I3.2[ for initially zero mean radial field all components 
of the mean magnetic field are ideally conserved. In Table [2] we show the com- 
ponents of the mean magnetic field, normalized to corresponding rms values. 
One can see that without proper mapping of the field fluxes the analytic con- 
straint is only fullfilled to roundoff error for B x . With the method described 
above one can also conserve the B z component to machine accuracy, while the 
error in the B y component can be reduced by about an order of magnitude. As 
can also be seen from Table [21 the maximum error in the solenoidal constraint 
is of the order 10~ 13 while the average error is as low as 5 x 10~ 15 . For the 
state of fully developed turbulence in the MRI-simulations (see section 14.31) 
the max./avg. values rise to 10~ 10 and 10~ 13 , respectively. 
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resolution 




24 3 


32 3 


48 3 


64 3 


/ R \ / Rrms 


a 


q snrio i fi 
y.ouue-iD 


l.UOZe-10 


1 /I9fia 1 
1.4fcZ0e-10 


O.Zooe-10 




D 


1 9Q9o 1 ^ 

i .zyze-io 


1.01Ze-lD 


9 nsfio i ^ 


O.DOoe-10 


/ R \ / Rrms 


a 




o.uooe-UD 


/i n9i o nR 


z.ouue-uo 




D 


Z.404ie-UD 


i ^97o nR 


/i nnrio n7 


1 nnso n7 
i.uuoe-u / 


/ R \ / Rrms 


a 


1 *?1 /la 


l.OUZe-UO 


a nfi 
O.0406-UD 


7 (\AF.a f\« 

1 .0406-UD 




D 


/I fiOSo 1 7 


A ^7/1 o 1 7 
41.0 / 4te-l / 


A fiHflo 1 7 


7 "^n/io 1 7 


mnv V7 . R /IRI 


a 




O.^^toe-lO 


7 /LzLfie 1 ^ 
1 .^^tue- io 


(J.U / ±e-±o 




b 


5.721e-13 


1.879e-13 


1.016e-12 


1.177e-12 


avg. V-B/|B| 


a 


5.493e-15 


5.469e-15 


4.575e-15 


3.964e-15 




b 


5.157e-15 


5.328e-15 


4.658e-15 


3.990e-15 



Table 2 

Convergence study of the relative error in the mean magnetic field components 
(normalized to rms- values) and solenoidal constraint: (a) without flux correction, 
(b) with mapped fluxes. 

In the following paragraph we analyze the epicyclic mode. As a real-life test 
case we have also performed simple MRI simulations with initially vertical 
magnetic field of zero net-flux. 



4.2 Conservation of epicyclic energy 

By performing two-dimensional shearingsheet simulations the error in the 
epicyclic energy is found to be independent of the mode amplitude. For con- 
stant excitation amplitude the error is growing linearly in time. From the 
last row of Table [3] one can see how the relative error per orbit decreases 
with resolution, reflecting the third-order convergence of the underlying time 
integration scheme. This is because due to the CFL stability condition, the 
timestep is linearly proportional to the spatial resolution. 

As already stated, our implementation of the source terms is integral part of 
the third-order integration scheme. For comparison we have also tested two 
conventional (operator split) methods for those terms: the first method (A) 
is very similar to the one found in the ZEUS-code and directly (i.e. forward 
Euler) integrates the Coriolis forces. The second method (B) treats the Coriolis 
term analytically in form of a rotation of the momentum vector. By expanding 
the trigonometric functions one can show that method A gives a first order 
approximation to method B. 
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32 64 
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128 



Fig. 2. Convergence of relative errors in epicyclic mode energy for method A (tri- 
angles), method B (squares), and the unsplit method (diamonds). Lines show least 
square fits with (logarithmic) slopes -1.79, -1.66, and -3.04. 



resolution 




24 2 


32 2 


48 2 


64 2 


96 2 


method A 


a 


2.776e-04 


1.637e-04 


7.943e-05 


4.877e-05 


2.535e-05 




b 


0.032180 


0.023989 


0.015896 


0.011886 


0.007900 


method B 


a 


7.387e-05 


4.497e-05 


2.301e-05 


1.473e-05 


8.135e-06 




b 


0.013782 


0.010275 


0.006810 


0.005093 


0.003385 


unsplit 


a 


4.045e-07 


1.695e-07 


4.862e-08 


2.023e-08 


5.345e-09 




b 


1.285e-09 


2.468e-10 


8.284e-09 


5.652e-09 


4.505e-09 




c 


4.038e-07 


1.695e-07 


4.874e-08 


1.747e-08 


2.587e-09 



Table 3 

Error growth-rates (a) and oscillation amplitudes (b) of the operator-split methods 
A (direct integration of Coriolis-forces) and B (exact rotation of momentum vectors). 
For the unsplit method we also show a simple linear fit (c). 



Both methods (in contrast to the unsplit one) lead to oscillations, with fre- 
quency 20, in the epicyclic energy. We fit the resulting curves with a function 
f(t) = at + bsin(2Qt). Table [3] shows the oscillation amplitudes and linear 
growth-rates in the relative error as a function of resolution. For reference we 
also include the fits for the unsplit scheme which show negligible oscillations. 
Figure [2] compares the rates of convergence for the described methods, clearly 
favoring the unsplit approach. 
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Fig. 3. Time evolution of volume averaged kinetic (thin lines) and magnetic (thick 
lines) energy density (top) and .R^-components of Reynolds- and Maxwell-stresses 
(bottom) for the non-conservative scheme (left) and the conservative scheme (right). 



4-3 MRI with vertical field of zero net- flux 



The most prominent application for the shearingbox model, of course, is the 
magneto- rotational instability (MRI). We chose model parameters according 
to previous simulations [E] with the old (non-conservative) NIRVANA-code, 
that is very similar to the widely used ZEUS-code. For simplicity we neglect 
stratification in this paper. A more sophisticated model, including stratifica- 
tion and radiative cooling has been implemented and will be employed in fu- 
ture work to explore MRI under conditions suitable to the interstellar medium. 

For now we use non-dimensional quantities, i.e., density is set to unity while 
the pressure is p = 0.5 x 1CT 6 , such that the sound speed over the box di- 
mension matches the angular velocity Q = 10~ 3 . The initial magnetic field is 
purely vertical and varies as B = B smijrx / L x ) z, resulting in a vanishing net 
vertical flux. The field amplitude Bq = 1.121 x 10~ 7 corresponds to a plasma 
parameter of (5 — 100 at the peaks of the sine-profile. 

We apply a box geometry of [—0.5,0.5] x [0,4] x [—2,2] with a resolution 
of 64x128x128 grid cells and perform computations on the new and the old 
code version to compare the conservative vs. the non-conservative scheme. The 
results are directly compared in Figure [3j The saturated stresses and energies 
are of comparable magnitude although the initial growth of the instability and 
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Fig. 4. Time evolution of volume averaged total (thick lines) and thermal (thin lines) 
energy for the non-conservative scheme (left) and the conservative scheme (right). 
For the former we also plot the work done by the boundary conditions (dashed line). 
All quantities are normalized by the initial gas pressure po- 

the breakdown of the channel-solution into chaotic turbulence differs quite a 
bit. The breakdown of the linear solution is related to parasitic instabilities 
that seem to be resolved differently in the two codes. The weaker initial peak 
for the conservative scheme is consistent with the results in |6J. Saturation 
amplitudes are compatible with earlier results [XlfT4"f6] . 

Although the codes compare quite well with respect to the kinematic quantities 
(see Figure [3]), there is a major difference concerning the thermalization of the 
extracted turbulent energy. As can be seen in the right panel of Figure H] the 
conservative code efficiently transfers the released energy into thermal energy, 
which leads to a rapid heating of the gas. The non-conservative code shows 
much less heating. Compared to the work exerted by the boundary-stresses^, 
that is shown as a dashed line (on top of the total energy) in the left panel of 
Figure HI there is a considerable amount of energy lost. 



5 Conclusions 



This paper described the implementation of a shearingbox environment for the 
flux-conservative/constraint-transport MHD code NIRVANA. To our knowl- 
edge this is the first time such an implementation has been developed for a 
central-type numerical scheme. We showed that shift-periodic boundary con- 
ditions accounting for a prescribed linear shear flow across the computational 
domain can be handled within a conservative framework. This comes about by 
a proper modification of the hydrodynamical numerical fluxes with the map- 
pings of y-momentum flux and total energy flux solely expressed in terms of 
known mapped fluxes and the velocity offset of the background shear flow. 
Such mappings are, in principle, independent of the underlying numerical 

4 This quantity is obtained by accumulating the rhs of Eq. (j27|) . 
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scheme, provided it is given in flux-conservation form. 

In contrast, our proposed mappings of the (face-centered) electric field fluxes 
contain explicit information on the magnetic field components at the boundary 
as well as on the characteristic speeds used by the numerical scheme. Never- 
theless, we demonstrated by numerical experiment that such a mapping is 
divergence-free, exactely conserves the x- and ^-component of mean magnetic 
field and the corresponding y-component up to truncation error but, obvi- 
ously, is scheme- dependent. An alternative method would be to directly map 
the edge-centered electric field fluxes, as those are the quantities which en- 
ter the constraint-transport equations. However, initial investigations of that 
approach have not been proven successful to date. 
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